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Abstract. Convection is the most important physical process that determines the structure 
of the envelopes of cool stars. It influences the surface radiation flux and the shape of ob- 
served spectral line profiles and is responsible for both generating and damping solar-like 
oscillations, among others. 3D numerical simulations of stellar surface convection have de- 
veloped into a powerful tool to model and analyse the physical mechanisms operating at the 
surface of cool stars. This review discusses the main principles of 3D stellar atmospheres 
used for such applications. The requirements from stellar structure and evolution theory to 
use them as boundary conditions are analysed as well as the capabilities of using helio- 
and asteroseismology to reduce modelling uncertainties and probing the consistency and 
accuracy of 3D stellar atmospheres as part of this process. Simulations for the solar surface 
made by different teams are compared and some issues concerning the uncertainties of this 
modelling approach are discussed. 
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1 . Introduction - 3D Stellar 
Atmospheres 

The inhomogeneous solar surface is charac- 
terised by large scale structures originating 
from the solar convection zone as well as from 
the complex solar magnetic field. Globally, 
the solar convection zone is dominated by an 
extreme density and temperature stratification 
with a contrast of 625,000 : 1 and 350 : 1, 
respectively, when comparing values for these 
quantities as found for the bottom of the so- 
lar convection zone with those found for its 
surface. The model data underlying these es- 
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timates are supported by helioseismology and 
solar surface observations (for further Refer- 
ences related to these estimates see Kimka 
2009). As a result the entire solar convec- 
tion zone comprises 20 pressure scale heights 
spread over 30% of the solar radius. In spite 
of the slow solar rotation rate the timescales 
on which convection transports mass and heat 
through the entire zone are hence comparable 
to those of rotation. This does not hold for the 
surface layers of the Sun. With average veloc- 
ities reaching 2 to 3 km s _1 the granules and 
downflow regions observed at the solar surface 
change over timescales of a few minutes to less 
than half an hour, two orders of magnitudes 
faster than the timescale of rotation on these 
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length scales. Similar holds for other types of 
stars with convection zones reaching the stel- 
lar surface provided that the surface gravity is 
high enough such that the pressure scale height 
at optical depths of order unity is just a small 
fraction of the stellar radius. 

These physical relations are the main back- 
ground for the "box-in-a-star" approach to 3D 
stellar atmospheres. In this case, a small vol- 
ume of the entire convection zone is consid- 
ered for which the conservation law equations 
for the densities of mass, momentum, and en- 
ergy, have to be solved. The hydrodynamical 
equations are coupled to the equation for radia- 
tive transfer, which have to be solved as well, 
since at the surface of a star the diffusion ap- 
proximation is no longer valid. Realistic hy- 
drodynamical simulations of such a system re- 
quire a tabulated equation of state and tabu- 
lated opacities computed for the chemical com- 
position assumed for the object to be studied. 
The simulations performed for this setting de- 
scribe the evolution of such a sample volume 
of the stellar atmosphere (and upper envelope) 
in space and time, as represented on a grid 
extending over the entire simulation box and 
computed for a finite number of time steps. The 
capabilities of this approach and its success- 
ful recovery of observational data such as sur- 
face intensity images in the visual or detailed 
spectral line profiles have been discussed in 
ma ny research papers and co nference reviews 
(cf. IStein & Nordlundlll998l as w ell as exam- 
ples^ and_2eferences inlXsplund 20071 ISteffenl 
20071, iK upka 2009, a nd the extensi ve review of 
Nordlund. Stein. & Asplundll2009l) . 

In the following we provide a short re- 
view on how 3D stellar atmospheres can be 
used for the modelling of stellar structure and 
for helio- and asteroseismology. We discuss 
requirements on model stellar atmospheres 
posed by these branches of stellar physics and 
some applications to illustrate the demands on 
the 3D models. More recently, 3D stellar at- 
mospheres have been used to compute physical 
quantities difficult to access by means of direct 
observations or even seismology, such as tur- 
bulent pressure or skewness of vertical veloc- 
ities below the visible solar or stellar surface. 
We provide some first results from a compar- 



ison of different 3D stellar atmospheres used 
to compute statistically averaged quantities re- 
quired in stellar structure modelling and seis- 
mology. We also discuss some issues related to 
the resolution of 3D models and conclude this 
paper with a summary. 

2. Stellar Modelling 

A solar structure model or a stellar evolution 
track depends on a variety of hydrodynamical 
processes including convection. Since the as- 
sociated hydrodynamical timescales are many 
orders of magnitudes shorter than the thermal 
and nuclear ones during most phases of stel- 
lar evolution, 3D hydrodynamical simulations 
of solar or stellar evolution are computation- 
ally too demanding. This restriction holds even 
for simulations of solar convection on just the 
hydrodynamical timescales because of extreme 
stratification (20 pressure scale heights from 
the top to the bottom of the zone). 

From a stellar structure modelling point 
of view convection models are required, 
among others, to compute stellar radii of stars 
which have convective envelopes. For this 
purpose tables of the mixing length parame- 
ter and turbulent pressure obtained from hy- 
drodynamical simulations of the surface con- 
vection zo nes would in principle be suffi- 
cient (see iLudwig. Frevtag. & Steffenl Il998l 
iTrampedach et al.ll 19991 and more recent work 
from the same authors for some examples). 
But for stellar evolution modelling such ta- 
bles are inadequate. In particular, these tables 
would not describe the nature of convection 
in the stellar interior, where we need to com- 
pute the amount overshooting of convective 
zones into radiative ones, the precise convec- 
tive efficiency in those transition regions, and 
other phenomena which occur in stellar interi- 
ors such as semi-convection and further vari- 
eties of double-diffusive convection. The cur- 
rently used analytical, semi-empirical formu- 
lae for the description of these physical pro- 
cesses suff er from inconsistencies (cf. the dis- 
cussions inlCanuto & Mazzitellilll99lHCanutol 
Il993t ICharbonnel & Zahnl 120071) . Hence, the 
development of improved calibration methods 
is highly rewarding as are more refined simu- 
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lations, which would be suitable for deducing 
more versatile convection models. 

Classical procedures for the computation 
of stellar models with convective envelopes 
such as calibrating the mixing length have 
some severe limitations. Among others, they 
are based on integral quantities (radius, lu- 
minosity, etc.). Unless secondary calibration 
methods are used (which in turn require stel- 
lar models), these quantities are known only 
for rather few, nearby stars with sufficient ac- 
curacy. In the end, such methods are an am- 
biguous probe of convection models. 

Methods based on detailed photospheric 
properties alone (such as comparing observed 
with computed spectral line profiles) require 
quantities such as the mean temperature to 
be properly modelled as a function of depth. 
They hence provide some extra information to 
reduce ambiguities in calibrations and mod- 
els tests, but they have only limited implica- 
tions on how convection is modelled in the 
stellar interior and can still remain inconclu- 
sive for photospheric co nve ction modelling (cf . 
iMontalban etall 120041 and iHeiter etall l2002l 
for example). 

3. Helioseismology and 
Asteroseismology 

Helio- and asteroseismology have opened a 
new window to look inside stars. They allow 
a much more thorough understanding of the 
mean structure and hydrodynamical properties 
of stellar interiors due to convection and rota- 
tion. Traditional "seismological tests" of con- 
vection models are based on frequency infor- 
mation obtained from helioseismological or as- 
teroseismological observations. But in spite of 
the constraints they provide on measurements 
of the depth of the solar convection zone, on 
the extent of deviations from a radiative tem- 
perature gradient in the layers right underneath 
the solar convection zone (overshooting), or on 
the chemical composition, these methods are 
still limited by ambiguities. One of the reasons 
is that the acoustic size of a resonant cavity 
can be the same for different model structures 
within the frequency domain which is accessi- 
ble to measurements with sufficient accuracy. 



This is the physical background for two con- 
tradicting explanations suggested for the dif- 
ferences between observed solar p-mode fre- 
quencies and predicted ones based on standard 
solar models dBasu & Antial [l995 found that a 
much steeper superadiabatic temperature gra- 
dient as obtained with the convecti on model by 
ICanuto & Mazzitellil I1991L [l992 reduces the 
discrepancy between observations and compu- 
tations, while a similar reduction was found 
when using averaged model structures from 3D 
hydro dynamical simulations of the surface lay- 
ers by iRosenthal et al] 1 19991 who pointed out 
the role of turbulent pressure and differences 
in the model structure originating from the 3D, 
inhomogeneous radiative transfer — s ee also 
ICanuto &~C hristensen-Dalsgaard 1998 for fur- 
ther discussion). 

Observational data on amplitude and time 
dependence (mode linewidth and lifetime) of 
the p-modes, i.e. the study of p-mode excita- 
tion and p-mode damping, provides the extra 
information necessary to remove such ambi- 
guities. Mode excitation occurs due to shear 
stresses and e ntropy fluctuations. Based on ear- 
lier work b y iGoldreich & Keelevl (I 19771) and 
Balmforth dl992|), different s emi-analytical 
models have been develo ped (s e elHoudek et all 
| 1999t ISamadi & GoupillboOll lHoudekll2002t 



Chapli n et al. I2005h . Through predictions of 



velocity and convective flux these models 
link observations with quantities required to 
construct stellar models. Models or numeri- 
cal simulations are used to compute eigen- 
functions and associated eigenfrequencies, the 
mean structure (density, entropy, . . . ), spatial 
and temporal correlations of velocity and en- 
tropy as a function of length scale, the filling 
factor (fraction of horizontal area covered by 
upflows), and the ratio of vertical to horizontal 
root mean square velocities. 

Observations (with some modelling) pro- 
vide quantities such as mode mass, the height 
above the photosphere where the modes are 
measured, the mode line width at half maxi- 
mum, and the mean square of the mode sur- 
face velocity. Sama di et alj (|2006) have used 
standard solar structure models to compute 
the model dependent input d ata f or the ap- 
proach of ISam adi & Goupil (2001) required 
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to calculate p-mode excitation rates. They in- 
vestigated the influence of non-grey atmo- 
spheres in comparison with the grey dif- 
fusion approximation and compared a stan- 
dard mixing l ength model of convection with 
the m odel of ICanuto . Goldm an. & Mazzitellil 
dl996l) . The standard model (grey, mixing 
length treatment of convection) underestimates 
excitation rates by up to an order of magni- 
tude. Convection treatment was found to be 
more important for the excitation rate cal- 
culation than line blanketing (non-grey at- 
mospheres), but eventually also the convec- 
tion m odel by ICanuto. Goldman. & Mazzitellil 
( 1996) combined with a non-grey atmosphere, 
as in Heiter et alj d2002l) . was found to yield 
a model structure that still underestimates p- 
mode excitation rates by up to a factor of 3. 

By comparison, numerical simulations 
were found to reproduce observed excita- 
tion rates under the assumptio n of a non- 
Gaus sian temporal correlation dSamadi et alj 
120031) . A new, semi-analytical model to com- 
put e p-mode excita t ion rat es was presented 
by IBelkacem et al.1 d2006al). It is based on 
the model of IGrvanik & Hartmannl d2002l) and 
iGrvanik eT al. (2005) for computing third and 
fourth order correlation functions in flows 
dominated by turbulent convection. This model 
accounts for asymmetry between up- and 
downflows (skewness of velocity and temper- 
ature). Due to mass conservation a flow with 
broad upflows (granules !) has higher veloci- 
ties in its downflows and since the input power 
for p-modes depends quadratically on the ve- 
locity, JheJoJaHnput power increases in this 
case. IBelkacem et al.l d2006bl) have shown that 
this model together with the assumption of a 
Lorentzian distribution function taken for the 
temporal correlations predicts solar p-mode 
excitation rates in agreement with observa- 
tions. We note at this point that some of the 
input for their model had to be taken from nu- 
merical simulations (for instance, the horizon- 
tal area fraction covered by upflow regions). 

The applicability of their new approach 
to stars other than the Sun was shown by 
ISamadi et alj d2008l) by computing mode exci- 
tation rates for a Cen A and comparing them 
with observational predictions. The model pre- 



dictions were found in agreement with the data, 
but it was also pointed out that a higher ac- 
curacy of asteroseismological measurements is 
still needed to provide sufficiently tight con- 
straints for the models, which is now possi- 
ble with the results from the CoRoT mission 
dAppourchaux et alJl2008t iMichel et al.ll2008l) 
and soon by the Kepler mission. This demon- 
strates that helio- and asteroseismology have 
the capability to falsify models and simulations 
where classical methods remain ambiguous. 

4. Comparing Simulations, 
Resolution Issues 

How well do the currently available numerical 
simulations of surface convection agree about 
the mean structure of the uppermost layers of 
stars and on the correlation functions required 
to interpret data from helio- and asteroseismol- 
ogy? To find an answer one should consider, 
among others, the following questions: to what 
extent is the flow influenced by the boundary 
conditions? How large is the influence of the 
domain size? How important is non-grey radia- 
tive transfer for the structure of the superadia- 
batic layer developed by convection zones that 
reach the stellar surface? Do different numer- 
ical methods and viscosity models change the 
large scale (coherently structured) part of the 
flow? 

In the following we discuss some first re- 
sults of a comparison of solar surface convec- 
tion simulations performed with four different 
simulation codes. Each of these 3D simulations 
is set up for the "box-in-a-star" scenario with a 
rectangular box covering only a small fraction 
of the stellar surface layers, with a constant 
surface gravity, and with periodic boundaries 
along the ho rizontal directions. Results for the 
CO 5 BOLD dFrevtag. Steffen. & Dorchl 120021) 
have kindly been provided by M. Steffen. This 
code uses a Roe-scheme for advection (i.e. 
a non-linear numerical viscosity) combined 
with a subgrid-scale viscosity dSmagorinskvl 
119631) . For the two simulations shown the 
vertical boundaries are open and the chem- 
ical composition assum ed is that one of 
lGrevesse&Sa"u^al dl998). The high resolution 
case shown assumes non-grey radiative trans- 
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fer (with 5 opacity bins), a box volume of 
11.2 2 x 3.1 Mm 3 , a grid of 400 2 * 165 points, 
a constant horizontal grid size of 28 km and 
a variable vertical one of 12... 28 km. The 
deep simulation assumes grey radiative trans- 
fer, a box volume of 11.2 2 x 5.2 Mm 3 , a grid 
of 200 2 x 250 points, a constant horizontal grid 
size of 56 km and a fixed vertical one of 2 1 km. 

Two further simulations have kindly been 
provided by FJ. Robinson and are based on 
the sim ulatio n code presented in [K im & Chan 
d!998l) and IRobinson etakl d2003l) (named 
CKS code in the figures). Its numerical scheme 
uses plain second order finite diffe rences with 
a subgrid-scale viscosity (Smago rinskvll 1963b 
where the coefficients are boosted near shock 
fronts. Both sets have closed vertical bound- 
aries and grey radiative transfer is assumed, 
with a 3D Eddington approximation. The sim- 
ulatio n with smaller dom ain size ("case D" 
from IRobinson et al.l 20031) assumes a chemi- 
cal composition as in Grevesse & Noelsl(ll993l) 
and has a box volume of 2.9 2 x 3 Mm 3 , a grid 
size of 58 2 x 170 points, a constant horizontal 
grid size of 50 km and a fixed vertical one of 
17.6 km. The "2008" sim ulation assumes the 
chemi cal composition of iGrevesse & Sauvall 
(119981) and has a box volume of5.4 2 x3.6Mm 3 , 
a grid of 1 17 2 x 190 points, a constant horizon- 
tal grid size of 46 km and a fixed vertical one 
of 19 km. 

Another simulation, performed with 
the code by IStein & Nordlundl (1 19981) . has 
kindly been provided by R. Samadi and 
K. Belkacem (it is essentia lly the one used 
by iBelkacem et al.l l2006al) . The advection 
scheme of this code is stabilized with hyper- 
viscosity and an extra viscosity near shock 
fronts. Open vertical boundary conditions are 
assumed and for this simulation the chemical 
co mposition is similar t o that one published 
by IGrevesse & Sauvall (Il998l) . A non-grey 
radiative transfer (with 4 bins) is assumed as 
well as a box volume of 6 2 X 3 Mm 3 , a grid 
of 150 3 points, a constant horizontal grid size 
of 40 km and a variable vertical grid with an 
average spacing of 20 km. 

Finally, we present results from a new sim- 
ulation performed wi t h the ANTARES code 
(iMuthsamet al.ll2007L |2009). Due to its high 



order advection scheme based on essentially 
non-oscillatory methods, i.e. a non-linear nu- 
merical viscosity, physical viscosity sources 
(radiative and molecular ones) are sufficient 
to stabilize the code for the photospheric lay- 
ers contained in the simulation. For the stellar 
interior their contributions are negligible and 
the numerical scheme is stable without adding 
subgrid-scale or artificial diffusivities. For this 
si mulation the chemic al composition presented 
in IGrevesse & Noelsl d 19931) is used together 
with a box volume of 6 2 x 3 Mm 3 , a grid of 
150 2 x 190 points, a constant horizontal grid 
size of 40 km and a constant vertical grid spac- 
ing of 16 km, and grey radiative transfer is as- 
sumed. 

In Fig. [1] we compare the horizontally av- 
eraged mean temperature T, plotted as a func- 
tion of the logarithm of the horizontally aver- 
aged mean density p, for all six simulation sets. 
The averaging was performed for each hori- 
zontal layer and each point in time followed 
by a time average. Since there is no inversion 
of the mean density as a function of depth in 
any of the simulations, there is a unique map- 
ping of the mean density to depth and radius 
for each simulation. In the interior all simu- 
lations agree very well with each other. Near 
the superadiabatic peak the maximum spread 
in mean density for a given mean temperature 
is +15% and is due to different metallicities, 
different resolution, and different treatment of 
radiative transfer schemes being used. For the 
lower photosphere all but one of the simula- 
tions agree again quite well (the exception is 
the "case D" model which does not extend far 
enough into the photosphere to produce real- 
istic temperatures in that region). For the mid 
and upper photosphere, which are not shown 
here, the simulations with grey radiative trans- 
fer yield higher average temperatures, although 
the detailed behaviour for the outermost layers 
also depends on the exact implementation of 
the boundary conditions. 

In Fig.|2]we compare the root mean square 
values of the horizontally averaged vertical 
velocity fluctuations of all six simulations. 
Despite the simulations with closed bound- 
ary conditions have w lms values which tend 
towards zero near the top, they behave very 
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Fig. 1. Horizontally averaged mean temperature T as a function of the logarithm of the horizontally av- 
eraged mean density p for the six simulation sets. The region shown includes the uppermost layers of the 
solar interior, the superadiabatic peak, and the lower photosphere. 



similarly to the simulations with open vertical 
boundaries for the layers at the solar surface 
up to the highest velocities, which are found 
around the superadiabatic peak. Only the shal- 
low "model D" computed with the CKS code 
and the lower resolution, deep model com- 
puted with the CO s BOLD code predict slightly 
lower maximum values. Further inside, both 
simulations with the CKS code predict higher 
values for w^, which drops rapidly near the 
closed bottom. In turn, the averages from the 
ANTARES code are very similar to those com- 
puted b y K. Belkacem and R. Samadi with the 
code bv lStein & Nordlundl(ll998l) from the so- 
lar surface close to the bottom of the simulation 
volume. Note the different behaviour of the 
simulation codes using open vertical bound- 
aries near the top of the simulation domain (in- 
creasing vs. constant Wrms). 

In Fig. [3] we compare the root mean 
square values of the horizontally averaged tem- 



perature fluctuations of all six simulations. 
Differences in the mid and upper photosphere 
originate from the non-grey radiative transfer 
assumed in some of the simulations. Except for 
the most shallow model, which has lower peak 
values, the root mean square fluctuations agree 
well among all of the simulations throughout 
most of the convection zone. Simulations with 
non-grey radiative transfer have larger peak 
values by only a few percent. Simulations with 
closed lower vertical boundary can be distin- 
guished by an increase of temperature fluctua- 
tions close to the bottom boundary layer. 

In Fig. |4] we compare the kurtosis K w of 
the horizontally averaged vertical velocity fluc- 
tuations for all six simulations. We note that 
again all simulations agree for the layers from 
the bottom of the photosphere till below the su- 
peradiabatic peak. However, further below the 
differences become large despite long-term av- 
erages performed over about 4 to 12 hours of 
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Fig. 2. Root mean square values of horizontally averaged vertical velocity fluctuations of all six simula- 
tions. 



solar time (about 15 to 50 convective turnover 
times) have been used. In particular, simula- 
tions with open vertical boundaries do not (!) 
agree more closely among each other when 
compared to simulations with closed vertical 
boundaries. This is an important caveat to re- 
member when using 3D stellar atmospheres to 
compute data for helio- and asteroseismology. 
Clearly, the higher order correlations are much 
more sensitive to the details of the implemen- 
tation of the boundary conditions than the vari- 
ables representing only the mean structure of 
the simulation domains. 

Since the simulations compared above 
have the same resolution within a factor of 
two, one might argue how they depend on 
numerical methods used and on resolution in 
general. None of the simulations compared 
above can be expected to resolve turbulence 
created by shear between up- and downflows, 
since in 3D one expects a horizontal and ver- 
tical resolution of less than 5 km to be re- 



quired for that dKupkal [2009). Questions re- 
lated to resolution include whether the spec- 
tral line profiles remain unaffected by turbu- 
lence generated through shear underneath the 
visible surface, what are the contributions of 
turbulence to Reynolds stresses and entropy 
for p-mode driving, or whether a better re- 
solved radiative cooling at the stellar surface 
does change the flow properties. Resolution ef- 
fects in 3D stellar atmosphere for grid spacings 
similar to the simulat ions discussed abo v e have 
been discussed by I Stein & Nordlundl d 19981) 
and iRobinson et alj (|2003l). a mong others. 



Recently, Muthsamet al. (2009) have demon- 
strated that for this resolution the diffusion of 
the advection scheme in the simulation clearly 
has an influence on the small scale structures at 
the visible surface. However, horizontal resolu- 
tions of 10 km with the least dissipative numer- 
ical scheme are necessary to finally reveal the 
highly turbulent nature of downdrafts under- 
neath the solar surface rather than just noticing 
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Fig. 3. Root mean square values of the horizontally averaged temperature fluctuations of all six simula- 
tions. 



indications for their turbulent behaviour. At the 
transition between laminar and turbulent flow 
the effective resolution of a numerical scheme 
is crucial, since only for even higher resolu- 
tion (i.e. in the turbulent regime) such numeri- 
cal details can be expected to no longer matter. 
The influence of shear driven turbulence on ob- 
servable quantities has yet to be investigated. 

5. Summary 

Hydrodynamical simulations have developed 
into a tool for the computation of realistic 3D 
stellar atmospheres of cool stars. Particularly 
for stars at or close to the main sequence they 
have successfully passed observational tests 
based on spectroscopy and traditional calibra- 
tion methods for convection models, but also 
tests based on helioseismology, which can dis- 
tinguish much more sensitively between dif- 
ferent physical models of the outer envelope 
of the Sun. Asteroseismology is now gradually 



gaining the same capability. Current numerical 
simulations of solar surface convection are ro- 
bust in their predictions of the mean solar sur- 
face structure (different boundary conditions, 
radiative transfer treatment, resolution) and 
pass observational tests which require mod- 
erate spatial resolution. However, a compari- 
son of higher order correlation functions for 
velocity (and temperature) demonstrates that 
these quantities are sensitive to the implemen- 
tation of boundary conditions. Interestingly, 
they could be tested by means of helio- and 
asteroseismology. Refined advection numerics 
and high spatial resolution show how turbu- 
lence is generated in down flows shrouded by 
the su rface layers (cf. also IStein & Nordlundl 
2000 for a first study of this question). 

Thus, 3D stellar atmospheres provide stel- 
lar evolution modelling with a tool to calibrate 
convection models and improve stellar evolu- 
tion calculations. Both spectroscopy and aster- 
oseismology will continue to help in defining 
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the region of applicability of numerical simu- 
lations of stellar convection. At the same time, 
questions related to mixing due to convection 
zones and their behaviour deep inside the stars 
remain open. A convection model taking these 
processes into account is still in demand. 
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